SkyView Subroutine

private subroutine SkyView(dem, azimuth, view)

Compute the maximum angle of sky obstruction along 16 directions

References:

Zhang, Y., Chang, X. & Liang, J. Comparison of different algorithms for calculating the shading effects of topography on solar irradiance in a mountainous area. Environ Earth Sci 76, 295 (2017). https://doi.org/10.1007/s12665-017-6618-5

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem
real(kind=float), intent(in) :: azimuth(:)
type(ViewingAngle), intent(inout) :: view(:,:)

Variables

Type Visibility Attributes Name Initial
logical, public :: check
real(kind=float), public :: delta
real(kind=float), public :: deltax
real(kind=float), public :: deltay
integer(kind=short), public :: i
integer(kind=short), public :: i1
integer(kind=short), public :: j
integer(kind=short), public :: j1
integer(kind=short), public :: l
real(kind=float), public :: length
real(kind=float), public :: topAngle
real(kind=float), public :: topAngleMax
type(Coordinate), public :: x0y0
type(Coordinate), public :: x1y1

Source Code

SUBROUTINE SkyView &
!
( dem, azimuth, view)
    
IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real)   , INTENT(in) :: dem
REAL (KIND = float), INTENT(in) :: azimuth (:)  !azimuth angles [rad]

!Arguments with intent(inout):
TYPE (ViewingAngle), INTENT(inout) :: view(:,:)

!local declarations
INTEGER (KIND = short) :: i,j,l,i1,j1
TYPE (Coordinate) :: x0y0, x1y1
REAL (KIND = float) :: delta, deltax, deltay
REAL (KIND = float) :: length
LOGICAL :: check
REAL (KIND = float) :: topAngle, topAngleMax

!----------------------------end of declarations-------------------------------
x0y0 % system = dem % grid_mapping
x1y1 % system = dem % grid_mapping
delta = dem % cellsize / 2.

DO i = 1, dem % idim
    DO j = 1, dem % jdim
        IF (dem % mat (i,j) /= dem % nodata) THEN
            DO l = 1, SIZE (azimuth)
                !get coordinate of starting point
                CALL GetXY(i, j, dem, x0y0 % easting, x0y0 % northing)
        
                deltax = 0.
                deltay = 0.
                topAngle = 0.
                topAngleMax = 0.
                DO 
                  !add increment in x and y
                  deltax = deltax + delta * SIN (azimuth (l) )
                  deltay = deltay + delta * COS (azimuth (l) )
                  x1y1 % easting = x0y0 % easting + deltax
                  x1y1 % northing = x0y0 % northing + deltay
                    
                  !retrieve local coordinate of new cell
                  CALL GetIJ (x1y1 % easting, x1y1 % northing, &
                              dem, i1, j1, check)
                  
                  IF (.NOT. check ) THEN
                      EXIT
                  ELSE IF (dem % mat (i1, j1) == dem % nodata ) THEN
                      EXIT
                  END IF
                  
                  !compute distance 
                  length = Distance (x1y1, x0y0)
                  
                  !compute topographic angle
                  topAngle = ATAN (  (dem % mat (i1, j1) - dem % mat (i,j) ) &
                                     / length )
                  IF (topAngle > topAngleMax) topAngleMax = topAngle
                END DO
                view(i,j) % angle (l) = topAngleMax
            END DO
        END IF
    END DO
END DO

RETURN
END SUBROUTINE SkyView